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ABSTRACT 

We present a direct numerical and analytical study of driven supersonic MHD 
turbulence that is believed to govern the dynamics of star-forming molecular 
clouds. We describe statistical properties of the turbulence by measuring the 
velocity difference structure functions up to the fifth order. In particular, the 
velocity power spectrum in the inertial range is found to be close to ~ k~^''^^, 
and the velocity difference scales as (|Am|) ~ The results agree well with 

the Kolmogorov-Burgers analytical model suggested for supersonic turbulence 
in [astro-ph/0 108300]. We then generalize the model to more realistic, fractal 
structure of molecular clouds, and show that depending on the fractal dimension 
of a given molecular cloud, the theoretical value for the velocity spectrum spans 
the interval [—1.74- ■ ■ — 1.89], while the corresponding window for the velocity 
difference scaling exponent is [0.42 ■■■ 0.78]. 

Subject headings: MHD: Turbulence — ISM: dynamics — stars: formation 



1. Introduction 



Observations and numerical simulations of gas motion in molecular clouds show rather 
complex distributions of velocity and density fields, indicating that the motion is turbulent. 
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Although the nature and the scales of the driving force of the turbulence can vary, it has been 
established long ago that this turbulence is highly supersonic, with Mach numbers varying 
from cloud to cloud and reaching up to of the order of 30 at scales ~ 100 pc (Larson 1979, 
1981, 1992; Falgarone & Phillips 1990; Falgaronc, Puget, & Perault 1992; Myers & Gammie 
1999; Padoan & Nordlund 1999, 2000; Ossenkopf & Mac Low 2000; Ostriker, Stone, & 
Gammie 2001). 

Star formation is, of course, ultimately due to gravitational collapse of small Jeans- 
unstable cores. However, as is seen in high-resolution numerical simulation, the initial density 
fragmentation, that leads to creation of such cores, may be mainly due to strong supersonic 
turbulence, i.e. may be explained to large extent without invoking the effects of self-gravity 
(Padoan et al. 2001). The process of star formation can therefore be divided into two stages 
that may be approached separately. In the first stage, the supersonic, turbulent motion of 
the interstellar gas develops shocks that interact with each other and with the turbulent 
flow, which results in the emergence of rather complicated density structures. This stage 
should be described statistically, in terms of probability distribution of density and velocity 
fluctuations, or in terms of their moments; attempts to understand the complicated picture 
of turbulence by studying each particular structure can hardly be successful. On the other 
hand, the second stage, dealing with collapsing cores, can be approached on the grounds of 
classical dynamics, and depends on the effects of gravity and of the specific environment (e.g., 
pressure, temperature, magnetic fields). These two stages are not independent, however, 
since the first one sets the distribution of initial conditions for the second one, for further 
discussion see, e.g., Padoan & Nordlund (1999); Klein, Fisher, & McKee (2000); Burkert 
(2001); Klessen (2001); Elmegreen (2001). Even more, with some suitable definition of 
density clumps, the first stage of density fragmentation leads to the distribution of clumps 
over masses that already resembles the observed initial mass distribution function of stars 
(Padoan & Nordlund 2000; Padoan et al. 2000). 

In this paper we address both numerically and analytically the first, "turbulent" stage 
of star-formation. So far, there has been no analytical theory predicting the statistical 
properties of supersonic interstellar turbulence, notwithstanding the fact that supersonic 
conditions have been inferred from observations for more than 20 years. In particular, the 
turbulence scahng relations, referred to as the Larson's laws, concern the scahng of velocity 
and density fiuctuations with respect to the size of the fiuctuations (Larson 1979, 1981). 
These scalings seem to vary for different clouds, but most observations suggest that the 
velocity difference scaling scatters around (|Aiip)^/^ ~ L°-^ (Falgarone, Puget, & Perault 
1992). The corresponding velocity power spectra are steeper than the Kolmog orov one Ef^ ~ 
which would correspond to (|AmP)^/^ ~ L^^^. As for the density fiuctuations, the 
scaling of the peak density of fiuctuations on scales L is close to p{L) ~ (Falgarone, 
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Puget, & Perault 1992). However, these results should be taken with a certain degree of 
precaution. The error bars of available observations arc rather large, and the systematic 
errors are sometimes unknown. Moreover, there are fundamental reasons that prevent one 
from restoring the three-dimensional velocity correlators from the measured two-dimensional 
projections (Ballesteros-Paredes & Mac Low 2001). 

The high-resolution numerical results that recently became available (e.g. Padoan, Nord- 
lund, & Jones 1997; Porter, Woodward, & Pouquet 1998; Mac Low 1999; Porter ct al. 1999; 
Stone, Ostriker, & Gammic 1998; Padoan & Nordlund 1999, 2000; Padoan et al. 2000) shed 
some light onto the physics of supersonic turbulence. In this paper, wc present numerical 
simulations of supersonic, super- Alfvcnic turbulence, driven on large scales in such a way 
that the sonic Mach number M is of the order 10 and the Alfvenic Mach number is of the 
order 3. An analysis of simulations with different Mach numbers will be presented elsewhere 
(Jimenez, Padoan & Nordlund 2001). We propose an analytical theory that explains the 
results of our numerical findings, and discuss its applicability to molecular clouds. Wc are 
interested mostly in the velocity correlators, although an application to the density statistics 
is presented as well. The next section analyzes the results of the numerical simulations. In 
particular, the velocity-difference structure functions are constructed. In section 3 we show 
that the observed features agree well with the recently proposed Kolmogorov-Burgers theory 
of supersonic turbulence. Section 4 discusses apphcation of the results to molecular clouds. 
Conclusions are presented in section 5. 



2. Scaling laws in numerical simulations of supersonic turbulence 

The numerical simulations were performed with 250'^ and 500'^ resolutions for MHD tur- 
bulence with an isothermal equation of state, using the same method and program as in 
Padoan, Nordlund, & Jones (1997), Padoan & Nordlund (1999), Padoan et al. (2000), and 
Padoan et al. (2001). The turbulence was driven on large scales by a solenoidal external force 
with 1 < A; < 2, where k — 1 corresponds to the size of the periodic box. The solenoidal 
character of the forcing is not crucial for the velocity scaling in the inertial region — this 
was checked by comparing with runs with mixed comprcssional and solenoidal driving. The 
solenoidal driving was chosen to provide a better 'boundary condition' for the inertial inter- 
val in k space, since it turns out that the comprcssional to solenoidal ratio tends to become 
small in the inertial range. The external force sustains the supersonic gas motion (M = 10) 
in the simulations. The motions would otherwise decay on a time scale of the order of a 
crossing time and become sub-sonic due to the dissipation in shocks (Stone, Ostriker, & 
Gammie 1998; Padoan & Nordlund 1999). The real forcing is probably due to a turbulent 
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cascade from large scales, driven by supernovae and superbubbles (Korpi et al. 1999; de 
Avillez 2000), 



The supersonic turbulence exhibits rather interesting properties that we summarize as 
follows. First of all, the spectra of both the potential, Uc, and the solenoidal, u^, components 
of the velocity field are steeper than the Kolmogorov spectrum k~^/^. These components are 
defined according to V ■ Ug = and V x Uc = 0. In Fig. 1 we plot the solenoidal spectrum 
weighted by A;^''^^, since the theory that we present below predicts for the velocity spectrum 
of turbulence Ek ~ k~^''^'^. Second, the divergence-free, solenoidal part of the velocity field. 
Us, is generated quite effectively by such turbulence, contrary to the two-dimensional case, 
where turbulence without pressure is mostly potential. This effect of vorticity generation is 
analogous to the magnetic dynamo effect existing in 3D and non-existing in 2D turbulence 
(also, the presence of a magnetic field can help generate vorticity, as discussed by Vazquez- 
Semadeni, Passot, & Pouquet 1996). We find that the compressible part accounts for only 10- 
20 percent of the intensity of the velocity field, see Fig. 1. The ratio 7 ~ (ul) / (u^) can thus be 
chosen as a small parameter of the turbulence. Third, the dissipative structures look like two- 
dimensional shocks rather than one-dimensional filaments or vortices as in incompressible 
turbulence. To demonstrate this, we plot in Fig. 2 a randomly chosen two-dimensional cross- 
section of the density distribution in a physical simulation domain. The filaments seen on 
the picture correspond to the two-dimensional shock structures. 

To quantitatively characterize the statistical properties of the turbulence we have mea- 
sured the so-called structure functions of the velocity field (Frisch 1995). These functions 
are defined as: 



where u is the component of the velocity field perpendicular or parallel to the vector L. 
According to the chosen component, the structure functions are called transversal or longitu- 
dinal, respectively. In the inertial interval, the structure functions obey scaling laws, and the 
exponents ({p) may be determined; it is usually expected that both transversal and longitu- 
dinal functions have the same scaling. The power spectrum of the velocity field is the Fourier 
transform of the second-order structure function, and may be expressed as Ek ~ k~^~'^^'^\ 

We performed measurements of the structure functions up to p = 5. Since the Reynolds 
number in our simulations was not large enough to observe good scaling behavior of the 
structure functions, the method of Extended Self-Similarity (ESS) was applied. Namely, 
instead of plotting the structure functions themselves, we plotted the ratios of their loga- 
rithmic slopes. As was discovered by Benzi, et al. (1993) and Camussi & Benzi (1996), such 
ratios exhibit rather good (and correct) scaling behavior, even in systems with moderate 




(1) 
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Reynolds numbers. Then, if one knows from the theory that the third order structure func- 
tion scales with C(3) = 1, or if any other scaling exponent ({n) is known from numerics 
or observations with good precision, one can obtain the scalings of all the other structure 
functions. To illustrate this procedure, in Fig. 3 we show the plots of the transversal struc- 
ture functions. Si . . . S5, their logarithmic slopes, and the ratios of the slopes vs that of S3. 
We chose to deal with the transversal rather than longitudinal functions since the flow is 
shear- dominated, and the transversal functions are therefore found to have a better scaling 
behavior. As illustrated by Fig. 3, the structure functions themselves do not exhibit good 
scaling behavior, but their relative scaling is well established in a rather large interval. As has 
been suggested by DubruUe (1994), the ratios of scaling exponents may be universal while 
the scaling exponents themselves may be not. The reason is that the physical length L may 
not be a good scaling variable, and the universal scaling of all structure functions holds not 
with respect to L, but with respect to a certain function ^{L) that depends on the Reynolds 
number and other parameters of the system. 

From Fig. 3 we conclude that the observed properties of the supersonic turbulence can 
hardly be predicted in the framework of any known model of strong turbulence. Indeed, 
the Kolmogorov model, or its generalizations to MHD, cannot be directly applied, since the 
turbulence has very small pressure and rather weak magnetic fleld. The model of turbulence 
without pressure, developed for the potential velocity field (the Burgers model), does not 
work either since the 3D turbulent flow generates vorticity. In the next section we present 
a theory of supersonic turbulence that explains the obtained spectra on the basis of the so- 
called She-Leveque model of strong turbulence. This explanation was suggested in (Boldyrev 
2001), and was motivated by recent successful application of the model to incompressible 
MHD turbulence (Biskamp & Miiller 2000; Miiller & Biskamp 2000). 



3. Kolmogorov— Burgers model for supersonic turbulence 

The analytical model hinges on the observation that in the inertial range the turbulence 
is mostly incompressible, obeying the Kolmogorov naive scaling laws of velocity fluctuations 
(see below) , while in the dissipative range it behaves as Biugcrs turbulence, developing shock 
singularities. A model that relates the dissipative structures to the velocity scaling in the 
inertial interval was developed by She & Leveque (1994) and She & Waymire (1995). As 
was pointed out by Dubrulle (1994), this model represents the velocity energy cascade as a 
log-Poisson process and can be obtained as a suitable limit of the so-called random (3 model 
of turbulence (see, e.g., Frisch 1995). In general, the She and Leveque model has three input 
parameters. Two of them are the exponents of the naive (i.e. non-intermittent) scalings of 



-6- 



the velocity fluctuation, ui I , and of the characteristic time of the energy transfer at this 
scale, ti l/ui l^. In our case these parameters are related, A = 1 — 9, although in general 
this does not need to be so, as e.g. in the Iroshnikov-Kraichnan model of incompressible 
MHD turbulence (Grauer, Krug, & Marliani 1994; Politano & Pouquet 1995). The other 
parameter is the dimension of the most singular dissipative structure, D. The She-Leveque 
formula, as generalized by DubruUe (1994), then reads: 



where H — 1 — A/{3 — D). Since in our case the velocity field in the inertial interval is mostly 
incompressible, one can imagine that the energy transfer is due to the Kolmogorov cascade, 
i.e., © = 1/3, A = 1 — © = 2/3, while the dissipative structures are shocks, not filaments, 
so D — 2. With such input parameters, the formula gives: 



For the first five structure functions the model gives: ({l)/({3) = 0.42, ^(2)/^(3) = 0.74, 
C(4)/C(3) = 1.21, ^(5)/C(3) = 1.40, in excellent agreement with the numerical results (cf. 
Fig. 3). To find the absolute scalings, we need to know the scahng of at least one of the 
structure functions. As we mentioned in Sec. 2, the Reynolds number is not large enough to 
see good scaling of the structure functions. However, it is known that the Fourier transform 
of the second order structure function, i.e., the velocity spectrum, can have larger scaling 
interval than the second order structure function itself. Thus, knowing from the measurement 
of the spectrum that ({2) ^ 0.7- ■ -0.8, we infer ({3) ~ 0.95 • • • 1.08, and the Kolmogorov 
relation, ^(3) = 1, may indeed hold in the inertial interval. 

Note that our formula (3) coincides with the formula derived in (Biskamp & Miiller 
2000; Miiller & Biskamp 2000) for incompressible MHD turbulence, where (^{3) = 1, and in 
this sense both systems belong to the same class of universality, in accord with the ideas 
put forward by Dubrulle (1994) and She & Waymire (1995). However, these two systems 
are completely different, and the Kolmogorov relation, ^(3) = 1, that is exactly proved 
in MHD, is not rigorously established in our case, but rather inferred from the numerical 
simulations. Strictly speaking, the scaling exponents can be different for these systems; it is 
their ratios that are reliably obtained from our numerical simulations and are believed to be 
universal. For instance, there is a possibility that these exponents can vary from cloud to 
cloud, depending on the equation of state, mechanisms of dissipation, Mach numbers, etc. 

For completeness, we would like to present here a cartoon model which leads to the 
She-Leveque formula. Our discussion is rather similar to the work of Dubrulle (1994). The 
reader not interested in the details of the derivation can safely skip to Sec 4. We start with 



C(p)/C(3) = © (1 - A) p + (3 - L>) (1 - E®^') , 



(2) 




(3) 
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the simple /3 model of turbulence, formulated in a way suitable for generalization of the 
results to the random (3 model (a detailed discussion of such models may be found in (Frisch 
1995)). Let us assume that our system is divided into boxes of size / and let us concentrate 
on the density of the velocity energy fiux over scales in each box, e^. Now let us divide each 
box into smaller boxes, with size IT, where F < 1. The number of boxes increases, but now 
we want to leave only a fraction (3 of all newly formed boxes, and since the energy fiux is 
constant, (e;) = const, the density of the energy flux in each box should be increased by 1/(3. 
After that, start decreasing the size of boxes again and repeat all the procedure exactly 
as in the previous step. One may introduce W — ej+i/ej, which takes the value 1/(3 with 
probability (3, and with probability 1 — (3 at each step, independently of previous steps. 
The fraction of space occupied by boxes decreases at each step, and one can easily check 
that eventually our boxes will cover the structure with fractal dimension 

D = 3-log(/3)/log(F). (4) 

This is the dimension of the structure to which the cascade converges. It is easy to see 
that (ef) - r^P\ where 

T{p)^\og{WP)/\ogr. (5) 

Let us now apply this model to incompressible turbulence. According to the Kolmogorov 
refined similarity hypothesis, ei scales as uf/l; we will denote this as ei ~ uf/l. So we have 

C(P) = P/3 + t(p/3) = p/3 + (3 - D)(l - p/3), (6) 

which gives, e.g., a velocity spectrum steeper than the Kolmogorov one. The (3 model 
produces the linear relation (6) which is not what is observed. 

We now generalize this model in the following way. First, let us assume that the steps 
of size changes are very small, i.e., F = 1 — x/{3 — D), and (3 — 1 — x, where x ^ 0. These 
expressions are chosen to satisfy (4) up to the first order in small x. To modify the previous 
model we now assume that we do not disregard any newly formed boxes, whose fraction 
would he 1 — (3 = X, but instead fill them with W = (3i, while the other boxes are filled 
with some factor W = (32- This procedure leads to a so-called random (3 model. Since we 
have to preserve (e;) = const, we have to impose the relation x(3i + (32{1 — x) = 1 which 
gives (32 = 1 + x{l — (3i). Using formula (5), we now derive r(p) = C{(3i — l)p + C(l — /3f), 
where we introduced the co-dimension C = 3 — D. If we now assume the Kolmogorov 
relation ei ~ uf/l, we immediately recover formula (2), where C(3) = 1 and (3i plays the 
role of E. This shows that E describes the degree of intermittency in the She-Leveque 
model. However, in the compressible case the Kolmogorov relation is not proven. We can, 
however, hope that if we change the length variable Z — > Z = ^(Z) in such a way that the 
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scaling exponent of the third-order structure function with respect to I becomes C(3) = 1, the 
scaUng of all the other structure functions will be restored as well. To achieve that we have 
to define the new scaling variable as I — ^{l) — {uf)/{ei). Also, instead of the Kolmogorov 
relation, wc assume e^/ (e) oc uf /{uf), which is suggested by the quadratic nonlinearity of the 
Navier-Stokes equation. Redefining the scaling r(p) as (ef) ~ r^P\ we now easily get: 

C(P)/C(3) = p/3 - C{1 - A)p/3 + C{1 - Pf/'), (7) 

which is the required generalization of the She-Leveque formula. The unknown parameters /3i 
and C can be obtained, e.g., from numerical simulations. In this formula it is not obvious 
that C is a co-dimension of the most singular dissipative structure, since the cascade was 
organized not with respect to physical length I, but with respect to /. 



4. Application to molecular clouds 

We now ask how adequate our numerical and analytical results are to real molecular 
clouds. The density structures of molecular clouds, as obtained from observations, are rather 
complicated, and one may question the validity of our simple model assuming shock singu- 
larities. In fact, observational relation of molecular cloud mass to the cloud size is close 
to M ~ L^, where D is estimated to be greater than 2 sec, e.g. Larson (1992), and also 
Elmegreen & Falgarone (1996), where the estimate D = 2.3 ± 0.3 was presented. If the 
density were concentrated in elongated shocks, the scaling power would be 2. A plausi- 
ble resolution lies in the assumption that over a range of scales (not available in numerical 
simulations due to limited resolution) dense structures are folded to form a complex distri- 
bution having fractal dimension D > 2; a complicated fractal (even multifractal) density 
distribution is indeed inferred from observations (Chappell & Scalo 2001). It seems natural 
to assume that the energy dissipation occurs mostly in shocks (Stone, Ostriker, & Gammie 
1998; Ostriker, Stone, & Gammie 2001), and therefore the dissipative singularities should 
have the same dimension. Formula (2) thus relates the fractal dimension of a given molecular 
cloud to scaling properties of the velocity field in this cloud. In principle, this dimension 
can be different for different clouds, depending on mechanisms of cooling, dissipation, struc- 
ture of magnetic fields, etc. Substitution of the fiducial value D — 2.3 in our formula (2) 
leads to even better agreement of the theory with recent observations, see Brunt & Heyer 
(2001), producing (|Au|) ~ L°-^^ and Ek ~ k~^-^^, where we assumed that the Kolmogorov 
relation C(3) = 1 holds. 

It is interesting to note that formula (2) also gives the upper boundary for the dimension 
of density distribution. The parameter E introduced in (2) measures the degree of intermit- 
tency of turbulence. For example, for E — > 1 one recovers the Kolmogorov, non- intermittent. 
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scaling ((p) — p/3. By its definition, < E < 1, and therefore, D < 2|. Taking into account 
that on small scales the dissipative structures take the form of two-dimensional shocks, one 
may argue that the lower boundary for the density dimension is -D = 2. We thus estab- 
lished a rather narrow window for D consistent with our theory. With the aid of formula (2) 
one then finds the corresponding intervals for the structure functions scaling exponents that 
are admissible in our theory. In particular, the first order exponent is expected to lie in 
the interval [0.42 • • • 0.78], and the second order exponent, related to the velocity spectrum, 
within [0.74 • • • 0.89]. It is interesting that the observed value, D — 2.3, is rather close to the 
upper boundary. In Figs. 4 and 5 we plotted the scahng exponents and C(2) as functions 
of D for 2 < < 2|. We see that the exponents are very sensitive to the value of D, and it 
is therefore very hard to infer ^'s from measurements of D this would require measuring D 
with high precision. More practical test of the theory would be plotting observed ^(2) versus 
C(l), our theory predicts that the values should lie on the line shown in Fig. 6. 



5. Conclusions 



The presented model proves to be rather successful in explaining the statistical prop- 
erties of supersonic turbulence. In particular, we were able to explain the observational 
scaling relations of the turbulence in molecular clouds. It would be highly desirable to 
obtain a clear physical picture underlying the Log-Poisson energy cascade assumed in the 
derivation of (2). Also, the relation of the density distribution to the statistics of the velocity 
field requires a better understanding. Such a relation may be derived analytically in a simpli- 
fied one-dimensional case (Boldyrev 1998; Boldyrev & Brandenburg 2001), although in the 
three-dimensional case the question is still open. However, in the derivation of formula (7) 
we explicitly constructed the multifractal distribution of the energy dissipation e, which can 
possibly be related to the multifractal distribution of the density, since the energy dissipates 
in shocks where the density is accumulated. And finally, in the present paper we did not 
try to address the problem of the initial mass distribution of stars (but see Padoan & Nord- 
lund 2000). The process of turbulent density fragmentation leads to creation of complicated 
density structures and sets the initial conditions for gravitationally unstable density clumps. 
The effects of self-gravity, and the dynamics of collapsing density cores fall beyond the scope 
of the present paper. All these are the questions for the future. 
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Fig. 1. — Power spectra of driven, super-sonic turbulence. The top panel shows the 
solenoidal power, compensated by A;^'^"^, and averaged over two turn-over times, after about 
three t.o. times from the start of the run, in a numerical experiment with resolution 250^ 
(plusses), using random driving at 1 < /c < 2. Additional data points (diamonds) are from 
a 500^ experiment, continued for about two tenths of a turn-over time from a 250^ snap 
shot, which is enough to estabhsh the extended large k solenoidal spectrum. The sohd hnes 
are least squares fits of the data in the range 2 < A; < 50 to functional forms a/{k^ -\- hk"^), 
where p is the power law index in the inertial range (1.75 for the 250^ case and 1.74 for the 
500^ case). Dashed lines show comparison slopes with spectral indices different by ±0.1, 
from 1.74. The bottom panel shows the average ratio of compressional to solenoidal power 
in the same experiments. 
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Fig. 2. — A random, two-dimensional cut through the physical simulation domain of a 
snapshot from a 250^ experiment. The cross-section of the density field shows filamentary 
structures that correspond to sheet-like shock density structures. 
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Fig. 3. — Transversal structure functions computed for p — 1,2,3,4,5 (correspondingly, 
from bottom to top in each panel). The first panel shows Log- Log plots of the structure 
functions. The second panel shows the differential slopes, C(l)) ■ ■ ■ > C(5)- The scaling range 
is very short, due to the limited Reynolds number. The third panel presents the ratios of the 
differential slopes to C(3), which exhibit excellent scalings, in agreement with the Extended 
Self- Similarity hypothesis. These ratios are well described by our formula (3). 



- 15 - 



0.75 - 




0.4 I I I I I I I I 

2 2.05 2.1 2.15 2.2 2.25 2.3 2.35 

Fig. 4. — Normalized scaling exponent of the first-order structure function, C(l)/C(3), as a 
function of the dimensionality, D, of the most singular dissipative structure, which can be 
close to the fractal dimensionality of the cloud. We use formula (2) with © = 1/3, A = 2/3 
and D changes from 2 to 2^. The results of our numerical simulations correspond to I? = 2. 



-16- 




2 2.05 2.1 2.15 2.2 2.25 2.3 2.35 

Fig. 5. — Normalized scaling exponent of the second-order structure function, C(2)/C(3), as 
a function of the dimensionality, of the most singular dissipative structure, which can be 
close to the fractal dimensionality of the cloud. We use formula (2) with © = 1/3, A = 2/3 
and D changes from 2 to 2^. The results of our numerical simulations correspond to I? = 2. 
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Fig. 6. — Normalized scaling exponent C(2)/C(3) (vertical axis) as a function of the nor- 
malized scaling exponent C(l)/C(3) (horizontal axis), as given by formula (2) with = 1/3 
and A = 2/3. The corresponding change of parameter D is in the interval [2; 2|]. The results 
of our numerical simulations correspond to I? = 2, i.e. C(l)/C(3) = 0.42 and C(2)/C(3) = 0.74. 



